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Abstract 

In this paper, we investigate the modification of our expressions developed for the 
modeMndependent data analysis procedure of the reconstruction of the (time-averaged) 
one-dimensional velocity distribution of Galactic Weakly Interacting Massive Particles 
(WIMPs) with a non-negligible experimental threshold energy. Our numerical simula¬ 
tions show that, for a minimal reconstructable velocity of as high as 0(200) km/s, our 
model-independent modification of the estimator for the normalization constant could pro¬ 
vide precise reconstructed velocity distribution points to match the true WIMP velocity 
distribution with a < 10% bias. 


1 Introduction 


Currently, direct Dark Matter detection experiments searching for Weakly Interacting Massive 
Particles (WIMPs) are one of the promising methods for understanding the nature of Dark 
Matter (DM) and identifying them among new particles produced at colliders as well as studying 
the (sub)structure of our Galactic halo [1, 2, 3, 4], 

In our earlier work [5], we developed model-independent methods for reconstructing the 
(moments of the) time-averaged one-dimensional velocity distribution of halo WIMPs by using 
the measured recoil energies directly. However, with a few hundreds or even thousands recorded 
WIMP events, only estimates of the reconstructed velocity distribution with pretty large statis¬ 
tical uncertainties at a few (< 10) points could be obtained. Hence, in order to provide more 
detailed information about the WIMP velocity distribution, we introduced the Bayesian analysis 
into our reconstruction procedure for concretely determining, e.g. the position of the peak of 
the one-dimensional velocity distribution function and the values of the characteristic Solar and 
Earth’s Galactic velocities [6]. 

In our Monte Garlo simulations presented in Refs. [5, 6], the minimal experimental cut-off 
energies of data sets to be analyzed are assumed to be negligible. For experiments with heavy 
target nuclei, e.g. Ge or Xe, and once WIMPs are heavy ( > 100 GeV), the systematic bias 
caused by this assumption should be acceptable. However, once WIMPs are light ( < 50 GeV) 
and a light target nucleus, e.g. Si or Ar, is used for reconstructing the WIMP velocity distribution 
fi{v), effects of a non-negligible threshold energy has to be considered and the estimate of the 
normalization constant of fi{v) would need to be modihed properly. Therefore, as a supplement 
of our earlier works, we consider in this paper the needed modihcation of the normalization 
constant of fi{v) for the general case with a non-negligible experimental threshold energy. 

The remainder of this paper is organized as follows. In Sec. 2, we hrst review the model- 
independent method for reconstructing the time-averaged one-dimensional velocity distribution 
of halo WIMPs by using data from direct DM detection experiments directly. Then, in Sec. 3, 
we develop the modification of the normalization constant of fi{v) for a non-zero minimal 
experimental cut-off energy step by step. Numerical results of the modihed reconstruction of 
the WIMP velocity distribution based on the Monte Garlo simulation will be given. A systematic 
bias in this model-independent modihcation will be discussed particularly. We conclude in Sec. 4. 
Some technical details for our analysis will be given in Appendix. 


2 Model—independent reconstruction of the one—dimensional 
WIMP velocity distribution 

In this section, we hrst review in brief the model-independent method for reconstructing the 
(time-averaged) one-dimensional WIMP velocity distribution by using experimental data, i.e. mea¬ 
sured recoil energies, directly from direct detection experiments. Detailed derivations and dis¬ 
cussions can be found in Ref. [5]. 


2.1 From the recoil spectrum 

The basic expression for the diherential event rate for elastic WIMP-nucleus scattering is given 
by [1]: 
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Here R is the direct detection event rate, i.e. the number of events per unit time and unit mass 
of detector material, Q is the energy deposited in the detector, F{Q) is the elastic nuclear form 
factor, fi{v) is the one-dimensional velocity distribution function of the WIMPs impinging on 
the detector, v is the absolute value of the WIMP velocity in the laboratory frame. The constant 
coefficient A is dehned as M = po0'o/2m^mr n; where po is the WIMP density near the Earth 
and (Jo is the total cross section ignoring the form factor suppression. The reduced mass rrii-^N is 
dehned by mr,N = + ^n), where is the WIMP mass and ttin that of the target 

nucleus. Finally, Umin is the minimal incoming velocity of incident WIMPs that can deposit the 
energy Q in the detector: Umin = ay/Q with the transformation constant 


a = 



( 2 ) 


and Vmax is the maximal WIMP velocity in the Earth’s reference frame, which is related to the 
escape velocity from our Galaxy at the position of the Solar system, Uesc- 

In our earlier work [5], it was found that, by using a time-averaged recoil spectrum dR/dQ 
and assuming that no directional information exists, the normalized one-dimensional velocity 
distribution function of incident WIMPs can be solved from Eq. (1) directly as 
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where the normalization condition: 
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has been used and thus the normalization constant A/” is given by 
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Here the integral goes over the entire physically allowed range of recoil energies: starting at 
Q = 0, and the upper limit of the integral has been written as cxd. Note that, the velocity dis¬ 
tribution function of halo WIMPs reconstructed by Eq. (3) is independent of the local WIMP 
density po as well as of the WIMP-nucleus cross section (Jq. However, not only the overall 
normalization constant A/” given in Eq. (5), but also the shape of the velocity distribution trans¬ 
formed through Q = jo? depends on the WIMP mass (involved in the coefficient a dehned 
in Eq. (2)). 


2.2 From experimental data directly 

In order to avoid some model dependence during giving a functional form for the recoil spectrum 
dR/dQ needed in Eqs. (3) and (5), expressions that allow to reconstruct f\{v) directly from data 
(i.e. measured recoil energies) have also been developed [5]. 

Consider experimental data described by 

Qn-^ <Qn,i<Qn + ^ : i = 1, 2, • • •, A^^, n = 1, 2, • • •, H. (6) 

Here the entire experimental possible energy range between the minimal and maximal cut-offs 
Qrain and Qmax has been divided into B bins with central points Qn and widths bn- In each bin, 
Nn events will be recorded. Since the recoil spectrum dR/dQ is expected to be approximately 



exponential, in order to approximate the spectrum in a rather wider range, instead of the 
conventional standard linear approximation, the following exponential ansatz for the measured 
recoil spectrum (before normalized by the exposure S) in the nth bin has been introduced [5]: 
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Here = N^/bn is the standard estimator for {dR/dQ)expt at Q = Qn, kn is the logarithmic 
slope of the recoil spectrum in the nth Q—bin, which can be computed numerically from the 
average value of the measured recoil energies in this bin; 


= (I) coth 

where 
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Then the shifted point Qs^n in the ansatz (7), at which the leading systematic error due to the 
ansatz is minimal [5], can be estimated by 
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Note that Qs,n differs from the central point of the nth bin, Qn- Finally, substituting the ansatz 
(7) into Eq. (3) and then letting Q = Qs,n, we can obtain that 


/'l,rec(ns,n) 


‘^Qs,n'^n 

— In F\Q) 

— kn 

[F^iQs,n)\ 

[dQ 

Q—Qs,n 


( 11 ) 


Here Vs,n = Oi^Qs,n, and the normalization constant M given in Eq. (5) can be estimated directly 
from the data by 
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where the sum runs over all events in the sample. 


2.3 Windowing the data set 


In order to reduce the statistical uncertainty on the velocity distribution reconstructed by 
Eq. (11) and some uncontrolled systematic errors caused by neglecting terms of higher pow¬ 
ers of Q — Qn, as well as to offer a reasonable number of reconstructable velocity points Vg^n 
of fi{v), it has been introduced in Ref. [5] that one can first collect experimental data in rel¬ 
atively small bins with linearly increased widths and then combining varying numbers of bins 
into overlapping “windows”. Thus, we set that the bin widths satisfy bn = bi + {n — 1)5, Hence, 
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Here the increment 6 satishes 6 = 2(^(5max — gmin — Bbi^/B{B — 1), B being the total number 
of bins, and Q (min,max) are the experimental minimal and maximal cut-off energies. Assume up 
to nw bins are collected into a window, with smaller windows at the borders of the range of Q. 



In order to distinguish the numbers of bins and windows, hereafter Latin indices n, m, ■ ■ ■ 
are used to label bins, and Greek indices /r, z^, • • • to label windows. For 1 < /i < nw, the /ith 
window simply consists of the first fi bins; for nw < Az < .B, the /ith window consists of bins 
/i — nw + 1, h — nw + 2, • • •, /r; and for i? < /i < i? + nw — 1, the /ith window consists of 
the last nw — {.l^ — B) bins. This can also be described by introducing the indices and n^+ 
which label the hrst and last bins contributing to the /ith window, with 




1 , 


/i — nw + 1; 


and 


for /i < nw, 
for /i > nw, 


= 


/i, 

B, 


for < B, 
for /i > 5. 


(14a) 


(14b) 


The total number of windows dehned through Eqs. (14a) and (14b) is evidently W = B + nw — ^, 
i.e. I < < B + nw — 1- 

For a “windowed” data set, one can easily calculate the number of events per window as 
N,= ^n, (15) 

n=n^-. 

as well as, the average value of the measured recoil energies 
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where is the central point of the /ith window. The exponential ansatz in Eq. (7) is now 
assumed to hold over an entire window. We can then estimate the prefactor as 
Wfj, being the width of the /ith window. The logarithmic slope of the recoil spectrum in the /ith 
window, kfj,, as well as the shifted point (from the central point of each “window”, Q^) can 
be calculated as in Eqs. (8) and (10) with “bin” quantities replaced by “window” quantities. 
Finally, the covariance matrix of the estimates of fi{v) at adjacent values of Vs,f_i = a^Qs,^ is 
given by^ 
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2.4 Numerical results 

In this section, we present reconstruction results of the one-dimensional WIMP velocity distri¬ 
bution before taking the modihcation of the normalization constant of fi{v). 

First of all, since the lighter the WIMP mass, the more problematic the non-negligible 
experimental threshold energy, a light input WIMP mass of = 25 GeV has been considered 

^Note that contributions involving the statistical error on the estimator for Af in Eq. (12) should in principle 
also, but do not be included here. 



in our simulations. As discussed and shown in Ref. [7], with 0(500) recorded events (in one 
data set), the (light) input WIMP masses can be reconstructed pretty precisely with only ~ 10% 
(a few GeV) statistical uncertainties (more simulation results can also be found in Ref. [6]). 
Thus, in our simulations demonstrated here and in the next section, the reconstructed WIMP 
mass involved in the coefficient a for estimating the reconstructed points as well as the 
normalization constant A/” has been assumed to be known precisely with a negligible uncertainty. 

As in Ref. [5], ’^®Ge has been chosen as our detector material for reconstructing /i(t;).^ As 
in Refs. [5, 6], the WIMP-nucleus cross section appearing in the expression (1) for the recoil 
spectrum dR/dQ has been assumed to be only spin-independent (SI), = 10“® pb, and the 

commonly used analytic form for the elastic nuclear form factor: 


fUQ) 


3ji(gRi) 

qRi 




(18) 


has been adopted. Here Q is the recoil energy transferred from the incident WIMP to the target 
nucleus, ji{x) is a spherical Bessel function, q = y/2m]<iQ is the transferred 3-momentum, for 
the effective nuclear radius we use Ri = with Ra — 1.2 fm and a nuclear skin 

thickness s ~ 1 fm. 

By taking into account the orbital motion of the Solar system around our Galaxy as well as 
that of the Earth around the Sun, the shifted Maxwellian velocity distribution of halo WIMPs 
has been given by [1, 5]: 
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Here vq — 220 km/s is the Solar orbital speed around the Galactic center, and Ve is the time- 
dependent Earth’s velocity in the Galactic frame [8, 1]: 
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with tp ~ June 2nd, the date on which the velocity of the Earth relative to the WIMP halo 
is maximal^. Additionally, a common maximal cut-off on the one-dimensional WIMP velocity 
distribution has been set as v^ax = 700 km/s. 

In our simulation shown in Fig. 1, the experimental threshold energy has been set as 
Qmin = 2 keV. Due to the maximal cut-off on the one-dimensional WIMP velocity distribu¬ 
tion, Vmaxi a kinematic maximal cut-off energy. 
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has to be considered. Since for our target '^®Ge it is only Qmax.kin.Ce = 52.65 keV, the maximal 
experimental cut-off energy has been set to be only Qmax = 50 keV. Meanwhile, since the lighter 
the WIMP mass, the steeper the expected recoil energy spectrum, &i = 5 keV width of the hrst 
energy bin has been used and the energy range between Qmm and Qmax has been divided into 
only four bins {B = 4); up to two bins have been combined to a window and thus four windows 

^Note that, while for a WIMP mass of 0(100) GeV, the transformation constant a defined in Eq. (2) is larger 
with ^®Si as the target nucleus than with "^^Ge, for a WIMP mass of 0(25) GeV, the transformation constant a 
with "^^Ge is larger. 

•^As usual, in all our simulations the time dependence of the Earth’s velocity in the Galactic frame, the second 
term of Ve(t), will be ignored, i.e. Ve = 1.05 uq is used. 



^®Ge, Q^in > 2 keV, < 50 keV, 500 events, = 25 GeV, 4 bins, up to 2 bins per window 



Figure 1: The reconstructed (rough) velocity distribution (black crosses) with a ^®Ge target for 
an input WIMP mass of = 25 GeV and an experimental threshold energy of Qmin = 2 keV. 
The vertical error bars show the square roots of the diagonal entries of the covariance matrix 
estimated by Eq. (17), whereas the horizontal bars indicate the sizes of the windows used for 
estimating fi,rec{vs,^) by Eq. (11). The solid red curve is the generating shifted Maxwellian 
velocity distribution with an input value of vq = 220 km/s. See the text for further details. 


(W = 4)^ will be reconstructed [6]. Additionally, we assumed that all experimental systematic 
uncertainties as well as the uncertainty on the measurement of the recoil energy could be ignored. 
5,000 experiments with 500 total events on average® in one experiment have been simulated. 

In Fig. 1, we show the reconstructed (rough) velocity distribution (black crosses) with a ^®Ge 
target. The vertical error bars show the square roots of the diagonal entries of the covariance 
matrix estimated by Eq. (17), whereas the horizontal bars indicate the sizes of the windows 
used for estimating /i,rec(G,/i) by Eq. (11). As a comparison, the solid red curve indicates the 
generating shifted Maxwellian velocity distribution with an input value of Vq = 220 km/s. It can 
be seen clearly that the reconstructed velocity distribution is around two times overestimated. 
This indicates in turn that, by using Eq. (12), the normalization constant W of the velocity 
distribution function fi{v) is underestimated as only the half! Remind that the experimental 
threshold energy used here is as low as only 2 keV and the corresponding minimal cut-off of the 
velocity distribution is 134.44 km/s (see also Fig. 2). 

^Note that the last window is neglected automatically in the AMIDAS code [9, 10], due to a very few expected 
event number in the last bin (window). 

^Note that, for our numerical simulations presented in this paper, the actual number of generated signals in 
each simulated experiment is Poisson-distributed around the expectation value. 

















3 Modification of the estimator for the normalization 
constant A/* 

In this section, we take into account the effect of a non-negligible experimental threshold energy 
(Qmin > 0) and introduce a model-independent modification of the estimator for the normaliza¬ 
tion constant Af step by step. 


3.1 Non—zero minimal cut-off velocity 

First, we consider the minimal cut-off of the velocity distribution due to the non-zero experi¬ 
mental threshold energy, nmin(Qmin) = 'i'min) i^se of the normalization condition (4). From 

Eq. (3), since v = by using integration by parts, we can obtain that 
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Here we define Q^ax = (Qinax, gmax.kin), the smaller one between the experimental and 
kinematic cut-off energies and can be understood as the upper bound of the recoil energy of the 
recorded events; Q(min,max) are the experimental minimal and maximal cut-off energies. Since the 
WIMP-nucleus scattering spectrum is expected to be exponential, the term of {dR/dQ)Q^Q^ 
has been ignored here, whereas 
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with ri = Ni/hi, an estimated value of the measured recoil spectrum {dR/dQ)expt {before the 
normalization by the exposure S) at Q = ^min- Finally, In{Qmin, ^max) can be estimated through 
the sum running over all events in the data set that satisfy Qa G [gmin, ^max]- 
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Once we neglect the (small) contributions from both v > Umax (< 0.5%) and v < and 
approximate the normalization condition by 
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a modihed estimator for the normalization constant can be given from Eq. (22) as 
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Remind that the hrst (extra) term in the bracket is caused by the non-zero minimal cut-off 
velocity and vanishes once the threshold energy is negligible (gmin — 0). 



^®Ge, Q^ip > 2 keV, < 50 keV, 500 events, = 25 GeV, 4 bins, up to 2 bins per window 



Figure 2: The reconstructed (rough) velocity distribution (black crosses) with the nor¬ 
malization constant estimated by Eq. (26). The vertical dashed blue line indicates 
f tnin (Qinin = 2 keV) = 136.44 km/s. All parameters are as in Fig. 1. See the text for detailed 
discussions. 

In Fig. 2, we show the reconstructed (rough) velocity distribution (black crosses) with 
the normalization constant estimated by Eq. (26). The vertical dashed blue line indicates 
V rr,\r, {Qru\u = 2 keV) = 136.44 km/s. Thus the meshed green area above = 136.44 km/s de¬ 
notes the integral on the left-hand side of the normalization condition (25), whereas the shaded 
light-red area blow has been neglected. 

Fig. 2 shows that the reconstructed velocity distribution with the modihed normalization 
constant given by Eq. (26) is strongly improved and could already match the true (input) 
velocity distribution (the solid red curve) pretty well. However, Fig. 2 shows also clearly that 
the contribution below the non-zero minimal cut-off velocity (the shaded light-red area) 
would still be non-negligible! 

3.2 Contribution below the non—zero minimal cut-off velocity 

The reconstructed velocity distribution shown in Fig. 2 matches almost perfectly to the true (in¬ 
put) distribution. However, as pointed out in the previous section, since not only the tiny contri¬ 
bution from the v > Umax area (< 0.5%) but also a pretty large part from v < = 136.44 km/s 

(the shaded light-red area) has been omitted, the normalization condition (25) with the integral 
over /i(u) only between and Umax could be considerably underestimated. Hence, an esti¬ 
mator for the area under the (true) velocity distribution function in the velocity range [0, 
should be given. 

In this section, we suggest a model-independent estimate for the area below As sketched 














in Fig. 2, for this aim, the value of the velocity distribution function aX v = has to be 
estimated and we approximate then the area of n < simply by a triangle. We start with 
the reconstructed recoil spectrum in the hrst Q—bin which can be given by 
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By using Eq. (3), an expression, similar to Eq. (11), for the value of the reconstructed velocity 
distribution function aX v = v*-,„ can be found as 
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Combining Eqs. (22) and (28), the integral over f\{v) between 0 and can be given by 
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Therefore, we can obtain the model-independent approximation for the normalization constant 
of the reconstructed WIMP velocity distribution as 
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Note that the modihed normalization constant M given here depends now on the estimates of 
ri and /ci.® Hence, the covariance matrix of the estimates fi,rec{vs,fj,) needs to be modihed as^: 
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®In fact, the normalization constant Af given in Eq. (26) depends also on the estimates of ri and ki. However, 
note that, we didn’t consider a modification of the covariance matrix used for estimating the statistical uncertainty 
bars shown in Fig. 2 as well as in the upper frames of Figs. 4 and Figs. 5, since Eq. (26) is only an intermediate 
product for our final expression (30). 

^Here we neglect again the relatively much smaller correlation between the uncertainties on Iq and ri, ki. 



^®Ge, Q^ip > 2 keV, < 50 keV, 500 events, = 25 GeV, 4 bins, up to 2 bins per window 



Figure 3: The reconstructed (rough) velocity distribution (black crosses) with the normalization 
constant estimated by Eq. (30) and the modified statistical uncertainties given by Eq. (31). The 
dotted magenta horizontal lines indicate the reconstructed velocity distributions including the 
correction of the theoretical estimate of the difference between the triangular approximation and 
the integral between 0 and (see the next section). All parameters are as in Fig. 1. 
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All derivatives of fi,rec{vs,^) to ri, ki and needed here are given in Appendix A.2. 

In Fig. 3, we show the reconstructed (rough) velocity distribution (black crosses) with the 
normalization constant estimated by Eq. (30) and the modified statistical uncertainties given by 
Eq. (31). It can be found that, due to the contribution of the first term in the bracket in Eq. (30), 
the reconstructed velocity distribution points are now a little bit underestimated. Meanwhile, 
the statistical uncertainties (vertical bars) given by Eq. (31) becomes now a bit larger, except 
the first one, which is significantly reduced. 

As a stricter check of our model-independent modification of the normalization constant A/", 
in Figs. 4, we increase the experimental threshold energy to Qmin = 5 keV. For a germanium 
detector, the corresponding minimal cut-off velocity of the incident WIMPs is now 215.72 km/s 
and thus the contribution from the area of u < u Tiin becomes much larger. Hence, the under¬ 
estimate of the normalization constant given by Eq. (26) and in turn the overestimate of the 
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^®Ge, Q^jn > 5 keV, < 50 keV, 500 events, rrij^ = 25 GeV, 4 bins, up to 2 bins per w/indow 



^®Ge, Q^jn > 5 keV, < 50 keV, 500 events, nij^ = 25 GeV, 4 bins, up to 2 bins per w/indow 



Figure 4: As in Figs. 2 (upper) and 3 (lower), except that the experimental threshold energy 
has been increased to Qmin = 5 keV. Note that the vertical scale of fi^redv) is different here. 
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Figure 5: As in Figs. 4, except that ^®Si is used as the target nucleus. Note that the vertical 
scale of fi^reciv) is different here. 






























reconstructed velocity distribution points are more clear. In contrast, by using Eq. (30) the 
reconstructed points can still match the true (input) velocity distribution pretty well. However, 
and more importantly, the systematic bias between the triangular approximation and the real 
value of the integral of the v < area becomes more obviously and problematic. We will 
discuss this drawback in more details in the next section. 

Finally, in Figs. 5, a lighter target nucleus ^®Si is used as the detector material. For a silicon 
detector with an experimental threshold energy of 5 keV, the corresponding minimal cut-off 
velocity of the incident WIMPs is 189.64 km/s. As the case with a '^®Ge target, the estimator 
(30) for the normalization constant M combined with the modihed covariance matrix given in 
Eq. (31) can not only correct the overestimated reconstructed velocity distribution very well, 
but also strongly reduce the large statistical uncertainty on the first reconstructed distribution 
point /i,rec(G,i), although the statistical uncertainties are (a bit) larger than those given with 
the Ge target (cf. Figs. 4). 


3.3 Bias of the estimator (30) for the normalization constant Af 

As revealed in Figs. 2 to 5, our model-independent triangular estimator for the integral over 
fi{v) between 0 and is somehow overestimated. Then the normalization constant M and 
the reconstructed velocity distribution points are in turn (a bit) underestimated. In this section, 
we discuss therefore the bias of the estimator (30) for the normalization constant M in more 
details. 

Taking the most commonly used shifted Maxwellian velocity distribution fi^sh{v) given in 
Eq. (19) as our theoretical assumption®, in Table 1 we give the estimated ratios of the triangular 
approximation (to the integral over fi{v) in the range between 0 and Umin) fo the integral itself: 





(32) 


as well as the fractions of the difference between the triangular approximation and the integral 
between 0 and (i.e. the little overestimated amount) to the integral in the entire velocity 
range between 0 and Umax: 


-/o"-^"/i(t^) dv 

Io^--fi{v)dv • ^ ^ 

20 different values of from 10 to 300 km/s for four most commonly used detector materials 
are given. 

As references, in Figs. 3 to 5 the reconstructed velocity distributions taking into account the 
corrections of the theoretical estimate of the difference between the triangular approximation 
and the integral over fi{v) between 0 and fii^) dv, have also been given as 

the dotted magenta horizontal lines. It can then be seen clearly that, with this (hnal) model- 
dependent correction for the normalization constant M, the reconstructed velocity distributions 
could match the true (input) distribution function very precisely: The tiny differences would 
totally be negligible, compared to the much larger statistical uncertainties given with 0(500) 
WIMP events. 

^Recently, several modifications of the Maxwellian velocity distribution have been introduced (see e.g. Refs. [11, 
12, 13]). However, as described in the papers, the significant differences between these (empirical) expressions 
and the shifted Maxwellian velocity distribution fi^sh{v) are only in the high-velocity tail. 



<in [km/s] 

Qmin [keV] 

* 

. V ■ 

A min 
^0 

Jr^^fi{v)dv 

28Si 

40 Ar 

76Ge 

i36Xe 


10 

0.014 

0.013 

0.011 

0.008 

1.4997 

0.0012 

20 

0.056 

0.054 

0.043 

0.031 

1.4987 

0.0094 

30 

0.125 

0.120 

0.097 

0.069 

1.4970 

0.0315 

40 

0.222 

0.214 

0.172 

0.123 

1.4946 

0.0742 

50 

0.348 

0.334 

0.269 

0.192 

1.4916 

0.1435 

60 

0.500 

0.482 

0.387 

0.276 

1.4877 

0.2452 

70 

0.681 

0.655 

0.526 

0.376 

1.4830 

0.3839 

80 

0.890 

0.856 

0.688 

0.491 

1.4775 

0.5635 

90 

1.126 

1.083 

0.870 

0.621 

1.4711 

0.7868 

100 

1.390 

1.338 

1.074 

0.767 

1.4637 

1.0551 

120 

2.002 

1.926 

1.547 

1.105 

1.4458 

1.7246 

140 

2.725 

2.622 

2.106 

1.504 

1.4234 

2.5495 

160 

3.560 

3.424 

2.750 

1.964 

1.3960 

3.4757 

180 

4.504 

4.334 

3.481 

2.486 

1.3634 

4.4152 

200 

5.561 

5.350 

4.298 

3.069 

1.3253 

5.2473 

220 

6.729 

6.474 

5.200 

3.713 

1.2816 

5.8242 

240 

8.008 

7.705 

6.189 

4.419 

1.2323 

5.9807 

260 

9.398 

9.042 

7.263 

5.186 

1.1777 

5.5473 

280 

10.900 

10.487 

8.423 

6.015 

1.1180 

4.3668 

300 

12.512 

12.039 

9.670 

6.905 

1.0537 

2.3101 


Table 1: The theoretically estimated ratios of the triangular approximation (to the integral 
over fi{v) in the range between 0 and ) to the integral itself as well as the fractions of the 
difference between the triangular approximation and the integral between 0 and (i.e. the 
little overestimated amount) to the integral in the entire velocity range between 0 and Umax- 
The most commonly used model for the Galactic WIMP velocity distribution function, /i,sh(^^) 
given in Eq. (19), has been used. 20 different values of from 10 to 300 km/s for four most 
commonly used detector materials are given here. The Qmin values have been estimated for the 
WIMP mass of = 25 GeV. 

























































































































































Figure 6: The theoretically estimated fraction of the difference between the triangular approxi¬ 
mation and the integral over /i(n) between 0 and (i.e. the little overestimated amount) to 
the integral in the entire velocity range between 0 and Umax as a function of 


Note however that, giving the correction of the systematic bias requires a theoretically pre¬ 
dicted velocity distribution function. For practical use without prior knowledge about the one¬ 
dimensional WIMP velocity distribution as well as for improving the simple triangular approxi¬ 
mation to the integral of the v < area (depending on the estimate of fi{v) at n = Umin); 
could use an iterative procedure with the Bayesian reconstructed velocity distribution function 
[6]. Nevertheless, considering the pretty large statistical uncertainties on the reconstructed veloc¬ 
ity distribution points as well as the much narrower statistical uncertainty band of the Bayesian 
reconstructed velocity distribution [6], the effect of ignoring the much smaller systematic bias 
with or even without the model-dependent theoretical corrections could not be significant (at 
least in the next few years). 

Besides of Table 1, in Fig. 6 we give the theoretically estimated fraction of the difference 
between the triangular approximation and the integral over fi{v) between 0 and (i.e. the 
little overestimated amount) to the integral in the entire velocity range between 0 and Umax as 
a function of n* for reader’s reference. 


4 Summary and conclusions 

In this paper, we investigated the modification of our expressions developed for the model- 
independent data analysis procedure of the reconstruction of the (time-averaged) one-dimensional 
velocity distribution of Galactic WIMPs with a non-negligible experimental threshold energy. 

It our earlier work [5], the experimental maximal and minimal cut-off energies have been 






















































































assumed to be large or small enough. Thus the sum over all recorded events in the data set can 
be used as the estimator for the integral over the one-dimensional WIMP velocity distribution 
function, which is in turn needed for the estimation of the normalization constant of the re¬ 
constructed velocity distribution. For experiments with heavy target nuclei, e.g. Ge or Xe, and 
once WIMPs are heavy ( > 100 GeV), the systematic bias caused by this assumption should be 
acceptable. However, once WIMPs are light ( < 50 GeV) and a light target nucleus, e.g. Si or 
Ar, is used for reconstructing the WIMP velocity distribution fi{v), effects of a non-negligible 
threshold energy has to be considered and the estimate of the normalization constant of fi{v) 
would need in turn to be modified properly. 

In this work, we derived at first the expression for estimating the integral over fi{v) between 
the minimal and maximal reconstructable velocities. Then we suggested the simple model- 
independent triangular approximation to the contribution below the minimal reconstructable 
velocity. Finally, by adopting the most commonly used shifted Maxiwellian velocity distribution 
function, the correction of the systematic bias caused by the use of the simple triangular ap¬ 
proximation has been given. Our numerical simulations presented in this paper show that, for a 
minimal reconstructable velocity of as high as 0(200) km/s, our model-independent modihca- 
tion of the estimator for the normalization constant (with or even without the model-dependent 
correction of the systematic bias) could provide precise reconstructed velocity distribution points 
to match the true WIMP velocity distribution with a < 10% bias. 

In summary, as a supplement of our earlier works on the (Bayesian) reconstruction of the 
WIMP velocity distribution function, we developed in this paper a model-independent modih- 
cation of the estimator for the normalization constant of fi{v) for the more general case with 
a non-negligible experimental threshold energy. This modihcation should not only be more 
suitable for our Bayesian reconstruction of the one-dimensional WIMP velocity distribution 
function [6], but hopefully also offer preciser information about Galactic Dark Matter for direct 
and indirect detection experiments and phenomenology. 
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A Formulae for estimating statistical uncertainties 

Here we list all formulae needed for the model-independent method for the reconstruction of 
the one-dimensional WIMP velocity distribution function described in Sec. 2 as well as for the 
modihcation of the normalization constant A/” by Eq. (30) and the modihed covariance matrix 
of the estimates of /i,rec('^s.n)- Detailed derivations and discussions can be found in Ref. [5]. 

A.l Formulae needed in Sec. 2 

First, by using the standard Gaussian error propagation, the expressions for the uncertainties 
on the standard estimator r„ and the logarithmic slope kn estimated by Eq. (8) can be given 



directly as 


cr r„ I = 


62 ’ 


and 




knkn/‘^ 


smh{knbn/2) 


-2 


a 


{Q - Qr, 


where 
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{Q-Qr. 


N„,-l 


{Q - QnY\n - Q - Q 
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n\n 


(Al) 


(A2) 


(A3) 


For replacing the “bin” quantities by “window” quantities, one needs the covariance matrix for 
Q — QaiIm, which follows directly from the dehnition (16); 

cov Qfilfj,-} Q Qi'liy'j 




N N 

-‘■'I1/ n=nv- 


Nn {Q\n ~ Q\i^ (Q\n ~ Q\i^ + A^n*^^ {Q ~ Qr. 


(A4) 


Note that, hrstly, /r < z/ has been assumed here and the covariance matrix is, of course, sym¬ 
metric. Secondly, the sum is understood to vanish if the two windows /i, v do not overlap, i.e. if 
n^+ < riu-. Moreover, similar to Eq. (Al), we can get 


cov(r^,r^) =- V Nn, 

n±;C- 

where ^Jt, < v has again been taken. And the mixed covariance matrix can be given by 


cov 


^l1 Q Qr 


1 «+ _ _ 

E Nn{Q\n-Ql 




(A5) 


(A6) 


Note here that this sub-matrix is not symmetric under the exchange of /r and u. In the dehnition 
of n_ and n+ we therefore have to distinguish two cases: 


Tl— 

Tl— ^/J,—5 ^1^+5 


if yU < z/; 
if yU > z/. 


(A7) 


As before, the sum in Eq. (A6) is understood to vanish if n_ > n+. 

Furthermore, the covariance matrices involving the estimators of the logarithmic slopes 
estimated by Eq. (8) with replacing n ^ fi, can be given from Eq. (A2) as 


cov {k^, kQ = klkl 


k^,b^/2 


smh{kf,bf,/2) 


-1 


kubul2 


sinh.{k^biy/2) 


and 


cov (r^, kQ = kl{l- 


X cov ^Q Qfj,\fj,: Q ) 

kybyl2 


sm\i{kvbu/2) 


-1 


cov 


I Q Qr 


(AS) 


(A9) 



A.2 Derivatives of the modified estimates /i,rec(^s,/x) 

The modified normalization constant M given by Eq. (30) depends on the estimates of ri and 
/ci, as the first reconstructed point of the velocity distribution given in Eq. (11), /i,rec('ys,i)- For 
modifying the covariance matrix of the estimates of fi{v), one needs thus to distinguish the 
H = 1 case from the other /i 7 ^ 1 cases. 

First, for the general 7 ^ 1 case, one has 


flMVs,,) = - 
a 


2 

a 


f (it* 1 -I- M T (n D* 

./l,recV^min/Vmin “T T JOWmin; Vr 

^AinF2(g) 

(Ikc^ 


n -1 


X 


[F^iQs,, 

(Qmin Qs,l) 


X 


F^iQmin) 

T -fo(f?min) Q 

2 g 


_d_ 

dQ 


Q — Qs,^j, 
2 / 


InF^(g) 


- h 


Q = Qn 


Qiain T 1 


-1 




Then the derivative of fi^redvs,^) to ri can be given as 
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(Alla) 


And the derivative of fl,rec{vs,^i) to ki is 
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Note that once gmin = 0, the second term in the bracket in the second line of Eq. (Alla) 
reduces to 1 and i9/i,rec('i's,/i)/'9ri as well as dfi^-recivs^Qj/dki become 0. Moreover, similar to the 
calculations done for the covariance matrix in Eq. (17), one can get 
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On the other hand, for the special /i = 1 case, we have 
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Then, similar to the calculations done in Eqs. (Alla) to (Al2b), it can be found that 
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